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Abstract. We investigate the thermodynamic properties of a dilute Bose gas 
in a correlated random potential using exact path integral Monte Carlo methods. 
The study is carried out in continuous space and disorder is produced in the 
simulations by a 3D speckle pattern with tunable intensity and correlation length. 
We calculate the shift of the supcrfluid transition temperature due to disorder and 
we highlight the role of quantum localization by comparing the critical chemical 
potential with the classical percolation threshold. The equation of state of the gas 
is determined in the regime of strong disorder, where superfluidity is suppressed 
and the normal phase exists down to very low temperatures. We find a T 2 
dependence of the energy in agreement with the expected behavior in the Bose 
glass phase. We also discuss the major role played by the disorder correlation 
length and we make contact with a Hartree-Fock mean-field approach that holds 
valid if the correlation length is very large. The density profiles are analyzed as a 
function of temperature and interaction strength. Effects of localization and the 
depletion of the order parameter are emphasized in the comparison between local 
condensate and total density. At very low temperature we find that the energy 
and the particle distribution of the gas are very well described by the T = 
Gross-Pitacvskii theory even in the regime of very strong disorder. 



PACS numbers: 05.30.Jp, 03.75.Hh, 03.75.Kk, 67.10.J 
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1. Introduction 

The dirty boson problem has become a central and fascinating subject in condensed 
matter physics starting from the first theoretical investigations more than 20 years 
ago [U [21 12] ■ The interplay between quantum degeneracy, interactions and quenched 
disorder in a bosonic system gives rise to a rich scenario that exhibits new and 
peculiar features compared to the much older problem of the metal/insulator transition 
with electrons [31 El E] ■ An important difference is that bosons can not rely on the 
Pauli pressure and repulsive interactions are crucial to avoid collapse in the lowest 
localized single-particle state. As a result perturbation schemes starting from the non- 
interacting model, that are most useful for fcrmions, are completely inappropriate in 
the case of bosons. 

Theoretical investigations, including quantum Monte Carlo simulations, have 
mainly addressed the problem of bosons on a lattice with on site bound disorder, 
the so-called disordered Bose- Hubbard model. In this case commensurability, i.e. the 
integer ratio of the number of particles to the number of lattice sites, plays a major role 
allowing for a superfluid/insulator (of the Mott type) transition also in the absence 
of disorder that is purely driven by interaction effects. Furthermore, depending 
on the value of the interaction strength, disorder can act in favor of superfluidity, 
by randomizing the insulating state close to the Mott transition, or in opposition 
to it by localizing almost free particles into single-particle levels. The insulating 
phases occurring in the two regimes of strong and weak interactions are respectively 
often referred to as the Bose glass, when interactions suppress superfluidity, and the 
Anderson glass, when interactions compete with disorder enhancing superfluidity [7] [8]. 
The disorder driven quantum phase transition occurring at T = has been investigated 
in a series of numerical studies both at incommensurate and commensurate densities 
and in various dimensionalities [7J GO [TU1 EJ H2J H31 1131 IS] • The picture emerging from 
these studies, together with the crucial role of interactions to stabilize the system, is 
that superfluidity is lost for strong enough disorder, leading to a gapless normal phase 
different from the incompressible Mott insulator. 

Random potentials in continuous systems have been considered using perturbative 
approaches based on the Bogoliubov theory [T5J [TCI El HH1 EI] • These methods are 
reliable when both interactions and disorder are weak and allow for the determination 
of the effect of disorder on the thermodynamic properties, including the fraction of 
atoms in the condensate, the superfluid density and other thermodynamic functions. 
Exact numerical methods have also been applied both at zero [21] and at finite 
temperature [5^1 In particular, the path- integral Monte Carlo simulations carried 
out at finite T addressed the problem of the elementary excitations [52] and of the 
transition temperature |23| of a Bose fluid in a random environment. In the case of 
the continuous-space liquid phase, disorder always acts against superfluidity, whereas 
interaction helps make the superfluid state more robust. For strong disorder and low 
temperatures one expects the system to enter an insulating phase (Bose glass) that 
is smoothly connected with the high-temperature normal phase existing when the 
disorder is weak. 

On the experimental side a large body of work was devoted to 4 Hc adsorbed 
in porous media, such as Vycor glass and aerogels. These studies investigated 
the behavior of the heat capacity and of the superfluid response (24] EH [26], as 
well as the dynamic structure factor [271 I28j as a function of temperature and 
filling. A suppression of the A transition is observed and the critical coverage for 
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the onset of superfluidity is determined as a function of temperature, however, no 
clear evidence is found of a compressible Bose glass phase. More recently the dirty 
boson problem has been addressed using ultracold atoms, which offer unprecedented 
control and tunability of the disorder parameters and of the interaction strength. 
Transport and phase-coherence properties of an interacting gas in disordered optical 
potentials are investigated and an insulating state is reached by increasing the 
strength of disorder j29l [30l [3TJ |32j [33j [3H [35] . A large experimental effort has also 
been devoted to the suppression of diffusion for non-interacting particles (Anderson 
localization) [31)1 |3"7] . 

In this work we report on a path-integral Monte Carlo (PIMC) study of an 
interacting Bose gas in the presence of correlated disorder produced by 3D optical 
speckles. This random potential is relevant for experiments and allows for an 
independent tuning of intensity and correlation length. By increasing the disorder 
strength, we find a sizable reduction of the superfluid transition temperature and the 
shift is larger for weaker interactions. We map out the normal to superfluid phase 
diagram, both in the chemical potential vs. disorder and in the density vs. disorder 
plane. For strong disorder and in the presence of small but finite interactions, the 
critical chemical potential varies linearly with the disorder intensity and is essentially 
independent of temperature and interaction strength, in agreement with the existence 
of a mobility edge separating localized from extended states. We also find that the 
critical chemical potential is much larger than the classical percolation threshold 
for the same random potential, implying that a major role is played by quantum 
localization effects. In the regime of strong disorder and for chemical potentials 
below the critical value, the equilibrium state is a highly degenerate normal gas. We 
investigate the thermodynamic properties of this phase, finding a T 2 dependence of 
the equation of state in agreement with the peculiar feature expected for the Bose 
glass phase. The effect of the disorder correlation length is discussed in detail and 
we show that a non-trivial behavior is obtained only when the correlation length 
is comparable with the mean interparticlc distance. At T = we also carry out 
calculations using the Gross-Pitaevskii (GP) equation and at finite T using a self- 
consistent mean-field approach based on Hartree-Fock theory and on the local density 
approximation. The results of the GP equation for the ground-state energy and the 
spatial distribution of particles are accurate even in the regime of strong disorder with 
short-range correlations. This conclusion might be useful in view of investigating the 
structural properties of the Bose glass phase. 

We consider a system of N identical bosons of mass m, subject to the random 
field Vdi s and interacting with a repulsive pairwise potential. The Hamiltonian of the 
gas takes then the form: 

N f h 2 \ 

i=l ^ ' i<j 

Interatomic interactions are modeled using a hard-sphere potential: 

V(r) = ( + n °° f r < a \ (2) 
w \ (r > a) , w 

where the diameter a coincides with the s-wave length of the two-body scattering 
problem. Furthermore, the system is in a cubic box of volume fl = L 3 with periodic 
boundary conditions. 
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Figure 1. (color online). Radial dependence (in units of the inverse momentum 
cut-off A) of the disorder spatial autocorrelation function I\ The solid (black) 
line refers to an average over many realizations of the random field, the (green) 
symbols correspond to a single realization. The (blue) line is a Gaussian fit. 



The structure of the paper is as follows. In section [2] we introduce the random 
potential and its statistical properties. In section[3]we discuss classical percolation for 
the speckle potential and we estimate the percolation threshold in 3D. Some details 
of the PIMC numerical method arc presented in section 01 In section Owe report our 
results on the superfluid transition: shift of the critical temperature, critical chemical 
potential and critical density. Most of these results were already presented in a 
previous publication of some of us |38j . In section|5]we introduce a mean-field approach 
based on the GP equation at T = and on a Hartree-Fock self-consistent theory at 
finite T and for long-correlated disorder. The low temperature thermodynamics is 
discussed in section[7J including the equation of state, the condensate and total density 
profiles and the behavior of superfluid density and condensate fraction as a function 
of temperature and interaction strength in the disordered phase. Finally, we draw our 
conclusions in section |8] 



2. Speckle potential 



The random external field we consider is the one produced by 3D optical speckles. 
The local intensity is obtained from the following expression [39] : 

2 



^dis(r) = V 



dk(p{k)W{k)e l 



where Vo is a positive constant and 
ip(k) = / dnp{r)e 



- ik-r 



(3) 



(4) 



is the Fourier transform of the complex field y(r), whose real and imaginary part are 
independent random variables sampled from a gaussian distribution with zero mean 
and unit variance. The function W(k) is a low-wavcvector filter defined as 

1 (k < 7rA) 

(k> ttA) . 



W(k) 



(5) 
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Figure 2. (color online). Typical shape of the speckle potential Vdisj with 
averaged value Vb = h 2 /ml 2 ., shown in the direction (0,0,1) of the simulation box. 




k B T/V 

Figure 3. (color online). Energy per particle of a classical non-interacting gas 
subject to the speckle potential. The homogeneous value E/N = 3kgT/2 is also 
shown as a reference. 



The random potential in equation §3§ is positive definite and the probability 
distribution of intensities is given by the normalized exponential law |40j 

P(V dis ) = L e -v^/v (6) 
Vo 

If the system's size is large enough the above defined disordered potential is self- 
averaging, i.e. spatial averages coincide with averages over different realizations. For 
a generic function f(V d i s ) of the disorder intensity one can thus write the following 
identity 

~ J drf[V dis (r)} = J q dV dis P(V dis )f(V dis ) = (f(V dis )) . (7) 

According to this property, the spatial average and the corresponding root-mean- 
square displacement of the speckle potential are both determined by the energy scale 



Dilute Bose gas with correlated disorder: A Path Integral Monte Carlo study 



6 



Vq- (Kiis) = Vq and yf (V^ is ) — (Vdis) 2 = Vq- The correlation length £ c of the random 
field is defined from the spatial autocorrelation function, 

r(r) = (V dis (r')V dis (r' + r)} - (V dis ) 2 (8) 

as the length scale for which Y(£ c /2) = T(0)/2. 

The above equations characterizing the speckle intensity field in 3D can be 
straightforwardly generalized to 2D and ID. In particular, in ID the autocorrelation 
function T(x) takes the simple form [50] 

/si^y 

\ irAx J 

and i c = 0.88/A. In 3D we find that the autocorrelation function Y(r) is well 
approximated by a gaussian (see figure [1]) and we obtain numerically £ c = 1.1/ A. 
It is important to notice that standard experimental realizations of optical speckles 
are 2D, i.e. the speckle pattern lies in the plane perpendicular to the propagation of 
the laser beam, featuring equal correlation lengths in the x and y planar directions 
and a much larger £ c in the z direction. We consider instead a 3D pattern, having 
the same correlation length in the three spatial directions. This random field can be 
realized, for example, by adding speckle patterns generated from different directions. 

The typical shape of the speckle potential Vdis is shown in figure [2j typical 
wells have size £ c and depth Vq . The energy h 2 /m£ 2 , associated with the correlation 
length £ c , and Vq provide the two relevant energy scales for the disorder potential. In 
particular, if Vq 3> h 2 /m£ 2 the random potential is classical in nature, with typical 
wells that are deep enough to sustain many single-particle bound states. The opposite 
regime, Vq <C h 2 /m£ 2 , corresponds instead to quantum disorder, where typical wells 
of size £ c do not have bound states and these can be supported only by rare wells of 
size much larger than £ c or with depth much larger than Vq. 

The root-mean-square intensity Vq and the correlation length £ c are the 
relevant parameters characterizing in general the various types of disorder. For 
example, the delta-correlated disorder, which has been considered in many theoretical 
investigations [TSJ [T8l [HI HQ] , is defined by the following autocorrelation function: 

<AF dis (r)A^ dis (r')> = nS(r - r') , (10) 

where A Vdis (r) = Vdis( r ) — (Vdis) is the displacement from the average. By 
approximating the speckle Y function ([SJ using a gaussian function, T(r) = 
V 2 e~ r l 2ry , with a = £ c / \/8 log 2 to recover the same half width at half maximum, one 
finds that the speckle field in the limit £ c — > reproduces a delta-correlated disorder 
with the strength k given by 

In our simulations the length scale £ c is typically ~ 100 times larger than the hard- 
sphere diameter a, allowing for a wide range of disorder intensities where interaction 
effects are well described by the s-wave scattering length and the details of the 
interatomic potential are irrelevant. The typical box size used in the simulations 
ranges from L ~ 20^ c to I ~ h0£ c . An indication of self-averaging of disorder for 
these values of L is provided by figure [T] where we show the comparison between the 
autocorrelation function Y averaged over many realizations of the random potential 
and the one corresponding to a single realization. 



Dilute Bose gas with correlated disorder: A Path Integral Monte Carlo study 



7 





Figure 4. (color online). Percolation in a 2D speckle potential: the left and right 
panel correspond respectively to an accessible volume (shown in red) below and 
above the percolation threshold. 



The self-averaging property (|7|) allows one to calculate the thermodynamics of a 
classical non-interacting gas. For example, the average energy per particle obtained 
from the spatial average of the disordered potential over the Boltzmann factor 

E_3_ fdrV dis (r)e- v ^/ k * T 

N~2 KB + Jdre-Wr)/fc B T > \ Li ) 
yields the following simple result 

§ = ^ + TTW- (13) 

In figure [3] we compare the above analytical prediction with the results obtained from 
a direct spatial integration using a typical size of the simulation box. The good 
agreement found shows that in our simulations the self-averaging property is well 
satisfied for non-trivial functions of the disorder intensity. 



3. Classical percolation 

In this Section we investigate the problem of the conducting/insulating transition in a 
speckle potential from the point of view of classical percolation. The relevant question 
is to determine the mobility edge of a classical particle subject to the random field [41] . 

Given the disordered potential Vdis (r) , the fraction of space accessible to particles 
of energy e is defined as the fractional volume where the reference energy exceeds the 
external field: 

m = i / dr. (14) 

The percolation threshold corresponds to the critical value 4> c of the fractional volume 
such that, if $(e) > $ c , there are infinitely extended volumes of allowed space and 
particles having energy e can move across the whole system. In terms of energy, the 
value e c determines the percolation threshold: 3>(e c ) = <S? C - It corresponds to the 
classical mobility edge separating localized states with energy e < e c from delocalized 



Dilute Bose gas with correlated disorder: A Path Integral Monte Carlo study 



8 



5 4 
10 L 



2JM0 3 5-10" 4 7.510" 4 10" 3 

<5 



Figure 5. (color online). Fraction of configurations of a 3D speckle pattern 
where percolation occurs in one of the spatial directions as a function of the 
accessible volume. Results are shown for two different system sizes. The estimated 
percolation threshold is 4(1) ■ 10 — 4 and is shown with the shaded area. 



ones with energy e > e c . In the case of speckles the function $(e) can be simply 
expressed in terms of the disorder intensity Vo using the property ([7]). One finds 

$(e) = 1 - e- £ / y ° . (15) 

The determination of the percolation threshold for lattice or continuum models 
is in general a difficult numerical task and the study of percolation has become a 
mature branch of statistical physics (see e.g. |H]). Our discussion here is limited to 
the calculation of <J> C and of the corresponding thercshold energy e c for the speckle 
potential. The role of dimensionality is crucial for this problem. In ID the unbound 
nature of the disordered field implies that, approaching the thermodynamic limit, 
potential barriers occur higher than any finite energy e and consequently e c = +00 
and <& c = 1. In 2D the percolation threshold of laser speckles was investigated 
experimentally [?3] using photolithography on a conducting film obtaining the value 
$ c = 0.407. This result was later confirmed by a numerical study [44] , To our 
knowledge there are no determination of $ c for 3D speckle patterns. 

We estimate the percolation threshold in the continuum by mapping the accessible 
and unacccssible regions of space on a finite grid. This is done by simply comparing the 
local value of the external field Vdi s (i"i) at the grid point with the reference energy 
e. One then investigates the percolation of accessible grid points of the corresponding 
matrix. In figure [4] we show two typical configurations in the case of 2D speckles: 
panel a) corresponds to $ = 0.30 well below the percolation threshold and panel b) 
to $ = 0.50 where percolation, i.e. existence of an uninterrupted path of accessible 
points across the whole system in at least one of the spatial directions, has clearly 
occurred. Two issues have to be considered with care: i) the size L in units of the 
correlation length £ c must be increased in order to extrapolate to the thermodynamic 
limit and ii) for a given size L the grid points must be dense enough. In figure [S] 
we plot the fraction of configurations exhibiting percolation in 3D speckle patterns 
as a function of the accessible volume $ and for different system sizes. An estimate 
of the threshold gives $ c ~ e c /Vb = 4(1) ■ 10~ 4 . For comparison, we estimated the 
percolation threshold in 2D obtaining the value $ c ~ 0.4 in agreement with previous 
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determinations [131 [31] . The value found for 4> c in 3D is rather small, suggesting the 
presence of many long valleys where particles with a tiny fraction of the energy Vo can 
freely move across the whole system. One should notice that a similar 3D continuum 
model, the percolation of voids between overlapping spheres (the so called "Swiss 
cheese" model), gives the much larger value $> c ~ 0.03 [45J [46] for the percolation 
threshold. 



4. PIMC method 



The partition function Z of a bosonic system with inverse temperature (3 = (fcgT) 1 
is defined as the trace over all states of the density matrix p = e~@ H properly 
symmetrized. The partition function satisfies the convolution equation 

p J ' p 

x J dRa... J dR M p(R,R 2 ,T)...p(K M ,PR,r) , 

where r = 0/M, R collectively denotes the position vectors R = (ri, r2, rjv), 
PR denotes the position vectors with permuted labels PR = (rpm,r P ( 2 \, ... ,r P r N \) 
and the sum extends over the N\ permutations of N particles. The calculation of 
the partition function in equation (|16[) can be mapped to a classical-like simulation 
of polymeric chains with a number of beads M equal to the number of terms of the 
convolution integral. In a PIMC calculation, one makes use of suitable approximations 
for the density matrix p(R, R', r) at the higher temperature 1/r in equation (|16[) and 
performs the multidimensional integration over R, R2,...,Rjvf as well as the sum over 
permutations P by Monte Carlo sampling [51]. The statistical expectation value of a 
given operator O(R), 

<°> = \jT\ E / dRO(R)p(R, PR, (3) , (17) 
• P J 

is calculated by generating stochastically a set of configurations {Ri}, sampled from 
a probability density proportional to the symmetrized density matrix, and then by 
averaging over the set of values {0(Ri)}. 

An approximation for the high temperature density matrix, which is particularly 
well suited for studies of dilute gases, is based on the pair-product ansatz [5T] 

^R', T )=n^,rs,r)n y{^ , ;i - m 

i=l i<j Prel\ T Vi l ii> T ) 

In the above equation p\ is the single-particle ideal-gas density matrix 

p 1 (r,r' J ,r) = (^-) 3/2 e-(-^W(-M, (19) 

and p re i is the two-body density matrix of the interacting system, which depends on 
the relative coordinates ry = — Yj and r J ■ = r ■ — r'j , divided by the corresponding 
ideal-gas term 

^(^ ) ^,r ) =( i ^-) 3/2 e-(--^)W(^). (20) 
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The two-body density matrix at the inverse temperature r, p re ;(r, r', r), can be 
calculated for a given spherical potential V(r) using the partial- wave decomposition 

Pre i(r,r',T) = — Y / (^+l)Pe{cos0) (21) 

POO 



where Pt{x) is the Legendre polynomial of order £ and 9 is the angle between r and 
r'. The functions Rk,t(r) are solutions of the radial Schrodingcr equation 

h 2 (d 2 R k j 2dR k j £{£+!) 



dr 2 r dr 



-Rk 



h 2 k 2 

V{r)R k ,i = R k ,i , (22) 



with the asymptotic behavior 

[2 sin(kr - en/2 + Sj) 
RkAr) = y " , (23) 

holding for distances r Ro, where i?o is the range of the potential. The phase shift 
St of the partial wave of order £ is determined from the solution of equation (|22|) for 
the given interatomic potential V(r). 

For the hard-sphere potential ([2]) a simple analytical approximation of the high- 
temperature two-body density matrix due to Cao and Berne |52| has been proven 
to be highly accurate [S3]. The Cao-Berne density matrix /0^ e f (r, r',r) is obtained 
using the large momentum expansion of the hard-sphere phase shift St — — ka + £tt/2 
and the large momentum expansion of the solutions of the Schodinger equation (|22|) 
Rk,t( r ) — \/2/7rsin[fc(r — a)]/r. This yields the result 

PSfO^V) = ! _ a(r + r')-a 2 

w a -[rr'+a 2 -a(r+7-')](l+cose)m/(2h 2 T) 

Simulations arc based on the worm algorithm |54j . which allows for an 
efficient sampling of permutation cycles. In this scheme one samples both diagonal 
configurations, contributing to averages of the type (|17p where all paths are closed, 
and off-diagonal configurations where one path is open. These latter configurations 
contribute to the one-body density matrix (OBDM) defined as 



"l 



= ~Y, I dr 2 ---dr N p(R,PR,(3), (25) 

D <* 



where rpm = r^. The long-range behavior of the OBDM determines the condensate 
density 

no = lim ni(r, r'). (26) 

r— r' | — >oo 

In our simulations the largest displacement of the OBDM we consider is |r — r'| = L/2. 
If the size L is large enough the number Nq of condensate particles can be written as 



/drm^r'), (27) 
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where r' is fixed by the constraint |r — r'| = L/2 and we perform an average over 
the solid angle. The quantity under the integral corresponds to the local condensate 
density at position r, which could be highly non uniform in the presence of a random 
potential. 

Beside the condensate density no, in the present study we consider also the 
supcrfluid density p s . The supcrfluid component is the part of the fluid that remains 
at rest when an infinitely slow movement is applied to the walls that contain the 
system. In the path-integral formalism, the supcrfluid fraction of a fluid contained in 
a box with periodic boundary conditions can be related [5T] to the fluctuations of the 
winding number via the equation 

Ps m(W 2 ) 



p 3h 2 f3N ' 
The winding number W is defined as: 



(28) 



N M 



W = ^ ^ (rj n+1 - rj n ) . (29) 

i—l m—1 

It is a topological property of the configuration. It counts the net number of paths 
that cross any plane perpendicular to one axis. The worm algorithm is particularly 
suitable to perform ergodic random walks that span all possible winding numbers 
since it extends the configurations space by including configurations with an open 
path. Only the Monte Carlo moves that modify the open path can efficiently change 
the winding number. 

We perform calculations both in the canonical (at fixed density n) and in the 
grand-canonical ensemble (at fixed chemical potential p) [54]. We supplement the 
worm algorithm with two additional Monte Carlo updates that change the particle 
number N. The first update adds one particle to the system by placing a closed path 
at a randomly selected position. The second update erases a randomly selected closed 
path. The acceptance probability of the first (second) update is fixed by the fugacity 
e 0n (by its inverse), by the change in the interaction energy due to the path to be 
inserted (erased) and by the factor (jfc?) ^at takes into account the density 

_ 3 

change and the normalization of the free particle propagator C = j3 / m) 2 . 

5. Superfluid transition 

The effect of disorder is to suppress both the superfluid and the condensate density. In 
figure [6] we illustrate the behavior of these two quantities when the disorder strength 
is increased. The value of p s /p and no/n in the absence of disorder is determined by 
the temperature T and by the value of the gas parameter no?. The figure shows a 
linear decrease of the supcrfluid and condensate components with increasing disorder 
in the classical regime Vq > h 2 /ml 2 .. An important question is to understand whether 
disorder also affects the superfluid transition temperature. The results are shown in 
figure [7j The transition temperature T c is expressed in units of 

-i 2/3 



2nh 2 
mkvt 



C(3/2) 



(30) 



the critical temperature of the non-interacting gas with C(3/2) — 2.612. In these 
calculations the value of the scattering length and of the disorder correlation length 
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°0 2 4 6 8 10 12 14 

V [h 2 /m^ ] 

Figure 6. (color online). Supcrfluid density and condensate fraction as a function 
of the disorder intensity Vb- The particle density is here no 3 = 1CT 4 and the 
temperature is T = 0.5T® . 




0.5- 

2 4 6 8 10 12 14 

V [h 2 lme 2 c ] 

Figure 7. (color online). Supcrfluid transition temperature as a function of the 
disorder strength for two values of the gas parameter no? . Open and solid symbols 
refer respectively to T c determined from the supcrfluid and from the condensate 
fraction. The dashed line is the prediction of Ref. I18I at no 3 = 10 -4 shifted by 
(T c — T®) in the absence of disorder. 



are kept fixed and for the latter we choose the value l c = 0.6n -1 / 3 . such that there 
is typically one particle in each sphere of radius l c : ~ 1. We show results 

corresponding to two different densities. The reported values of T c in the absence of 
disorder are taken from Ref. [55] . At the density corresponding to the gas parameter 
na 3 = 10~ 4 we find no appreciable change of the transition temperature compared 
to clean systems by increasing the disorder strength up to Vb ~ h? jml\. For larger 
intensities we find a sizable shift that is well described again by a linear dependence 
in Vo- For a given strength Vq the reduction of the transition temperature is enhanced 
for smaller values of the gas parameter, consistently with the instability of the ideal 
Bose gas in the presence of disorder |16j . The value of T c is extracted from the results 



Dilute Bose gas with correlated disorder: A Path Integral Monte Carlo study 



13 





Figure 8. (color online). Scaling behavior of the superfiuid density for different 
system sizes and different realizations of disorder. The value of the gas parameter 
is na 3 = 10~ 4 . 
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Figure 9. (color online). Critical chemical potential (shifted by Vo) as a function 
of the disorder strength for different values of the temperature (in units of h 2 /m£ 2 ) 
and of the scattering length. The grey shaded area denotes the superfiuid phase 
and the (pink) solid line corresponds to the classical percolation threshold. 



of the superfiuid fraction p s j p (p = mn is the total mass density), corresponding to 
systems with different particle number N, using the scaling ansatz 

N^ 3 Ps (t, N)/p = /(iiV 1 / 3 ") = /(0) + /'(0)t7V 1/3iy + ... . (31) 

Here, t = (T — T c )/T c is the reduced temperature, v is the critical exponent of the 
correlation length ~ t~ v ', and f(x) is a universal analytic function, which allows 
for a linear expansion around x — 0. The validity of the scaling behavior (|31[) is shown 
in figure [51 where the effect of different realizations of the random potential is also 
shown. The quantity iV( 1+I? '/ 3 no/n, involving the condensate fraction n^jn and the 
correlation function critical exponent r\ = 0.038 of the XY-model universality class, 
is also expected to obey a scaling relation of the form (|3lj) . For all reported disorder 
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Figure 10. (color online). Spatial dependence of the OBDM for two values of 
the chemical potential slightly below and above p, c . Here kgT = 0.13S 2 /m£ 2 and 
a/£ c = 0.016. Two different system sizes are used to check the role of finite-size 
effects. 
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Figure 11. (color online). Critical density as a function of the disorder strength 
for different values of the temperature (in units of h 2 /ml^,) and of the scattering 
length. The critical density separates the supcrfluid phase (above n c ) from the 
normal phase (below n c ). The horizontal arrows indicate the critical value n° of 
the non-interacting gas. The (blue) triangle corresponds to the same temperature 
but a larger scattering length compared to the (green) squares. Inset: Density 
dependence of ps/p [(pink) squares] and no/n [(blue) circles] for the values 
Vq = and V = 6Ah 2 /me' 2 of the disorder strength. Here k B T = 0.13ft 2 /"^ 
and a/£ c = 0.016. The vertical arrow indicates the corresponding value of the 
degenerate density n^. 
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strengths Vo, the extracted value of the critical exponent v is compatible with the 
result v = 0.67 corresponding to clean systems [55]. It is worth noting that the values 
of T c , obtained from the scaling law of the supcrfluid p s / 1 p and of the condensate 
fraction no/n, coincide within our statistical uncertainty (see figure [7]). The presence 
of disorder reduces the quantum derealization of particles occupying the deepest 
wells of the potential and, consequently, their contribution to the supcrfluid behavior. 
Superfluidity takes place when the degeneracy condition is met for the effectively 
smaller density of "dclocalized" particles, resulting in a suppressed value of T c . In 
Ref. [T5] the shift 5T C = T c — of the superfluid transition temperature is calculated 
using a perturbative approach for the ^-correlated disorder (AVdis(r)AVdis(r')) = 
nS(r — r'), where AVdis(r) = Vdi S (r) — (Vdis)- The T c shift is found to be quadratic in 
K, implying for our speckle potential that ST C /T° = (m 2 Vg tl / ^Jwtf) 2 / [2(12 log2) 3 ], 
where we used a gaussian fit to the radial dependence of the autocorrelation function 
r and considered the limit l c — > 0. We report this prediction in figure [7J (we also 
add the interaction contribution not accounted for by Ref. [IB], so that in the clean 
case an exact result is reproduced). Our data in the regime of very weak disorder 
do not have enough precision to allow for a quantitative comparison and diverge 
from the theory before 5T C /T® becomes appreciable. The effect of disorder on the 
critical temperature of a hard-sphere gas was also investigated using PIMC methods 
in Ref. [23] where, however, no significant reduction of T c was reported. For stronger 
intensities of disorder, the calculation of T c becomes increasingly difficult, since the 
dependence on the realization gets more important (see figure [8]) and larger systems 
are needed in order to have a satisfactory self-averaging of the random potential. 

In figure [9] we report results for the critical chemical potential p c obtained from 
calculations carried out in the grand-canonical ensemble. A small change of p around 
p c translates into a drastic change in the long-range behavior of the OBDM (see 
figure [T0| : for \x < p c the OBDM decays to zero and corresponds to a normal 
phase, for p, > p c the OBDM reaches a constant value characteristic of the supcrfluid 
state. If interactions are small but finite, we also find that the value of p c is 
essentially insensitive to a change of temperature and of interaction strength. For 
weak disorder, this result is accompanied by a very small critical density (see figure ITT]) 
and corresponds to a renormalization of p c due to disorder in an extremely dilute 
gas. For strong disorder, it is instead consistent with the picture of a mobility 
edge, separating localized single-particle states from extended ones, that depends 
only on the parameters of the random potential. In this latter regime we find a 
linear dependence of p c as a function of Vo, in agreement with the qualitative T = 
prediction of Rcfs. [57] [56] in the case of classical disorder. The figure also shows 
the classical percolation threshold p = e c , whose value for the speckle potential has 
been determined in section [3J One should notice that in the whole range of disorder 
intensities the critical chemical potential is significantly larger than e c as a consequence 
of quantum localization effects. In fact, in terms of a mobility edge, classical particles 
with energy larger than e c can freely move across the entire system, while in the 
quantum world extended states appear only for significantly larger energies bound by 
the inequality e > p c . 

To conclude the study of the critical behavior, we analyze the dependence of 
the critical density n c on the intensity of the random potential. The calculations are 
carried out in the canonical ensemble at fixed temperature and scattering length. The 
method used to determine n c is shown in the inset of figure 111! For a given value 
of Vo one increases the density and calculates the superfluid p s /p and the condensate 
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Figure 12. (color online). Chemical potential at T = from the GP equation as 
a function of the disorder correlation length. The value of the gas parameter is 
na 3 = 10 _ 6 and the disorder strength is Vo = 10k gT^. The limits of small and 
large £ c are shown with horizontal lines. 



fraction no/n. The results are then fitted by a power-law dependence p s j p ~ (n — n c ) u 
and no/n ~ (n — n c ) v ^ 1+ri ' for n > n c , where the proportionality coefficients are 
expected to be non-universal parameters. In the inset of figure [TT] we show the results 
corresponding to a configuration without disorder (Vo = 0) and with strong disorder 
(Vo = 6 Ah 2 1 'm£ 2 ) . The reported values are averaged over a few realizations of the 
random potential and their scatter gives an idea of the relevance of this effect. For the 
small value of the scattering length used here, the critical density at Vq = coincides 
with the non-interacting result n° c = £(3/2)(mfcBT/27r/i 2 ) 3 / 2 , while for the large Vb 
one finds that n c is about a factor of eight greater than n°. More comprehensive 
results are shown in figure [11] where n c is estimated from the supcrfluid fraction, 
which is less sensitive to finite-size effects. The results clearly show an increase of the 
critical density as a function of Vq, from the non-interacting degenerate density riP c 
up to values ~ 20 times larger. It is also worth noticing that for strong disorder, an 
increase of the scattering length a is accompanied by a decrease of n c resulting in a 
constant value of the critical chemical potential (see figure . 

6. Mean-field approach 

A simple description of the thermodynamic properties of disordered systems can be 
provided in terms of a mean-field approach. At T = this approach is based on the 
solution of the Gross-Pitacvskii (GP) equation for the order parameter in the random 
external field and it yields quantitatively reliable results for both the chemical potential 
and the density profiles. At finite temperature the mean-field theory can be efficiently 
applied in the case of random potentials with exceedingly long-range correlations where 
the local density approximation holds valid. 
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Figure 13. (color online). Density dependence of the chemical potential at T = 
from the GP equation: (red) solid line. The (green) dashed line corresponds 
to the Thomas-Fermi limit fj, = \/2gnVo and the (blue) dotted line to the 
opposite limit of short correlation length fj, = gn + Vq. The disorder strength 
is Vq = 3.3 X I0~^h? /ma 2 and at the density ria 3 = 10~ 6 corresponds to the 
value used in figur dl2l 



6.1. Zero temperature 



The relevant equation is the stationary GP equation for the order parameter ip(r) in 
the presence of the random potential 

h 2 



2m 



V" + Vdis(r)+ ff |V(r)|' 



V>(r) = (j,ip(r) 



(32) 



where g = 4 7r ft 2 a/m is the coupling constant. Within the GP approach one does not 
distinguish between order parameter and particle density and the following identity 
holds valid: n(r) = |?/;(r)| 2 - The GP equation can be obtained from the following 
energy functional [37] 



E[ip] 



d 3 rV>*(r) 



-|-V 2 + F dii 
2 m 



V>(r) 



+ 2<#(r)| 4 



using the variational ansatz 



6 

Si? 



d 3 r\ij)\ 



= 



(33) 



(34) 



that corresponds to finding configurations that minimize the energy functional (j34J) 
with the normalization constraint J d 3 r\tp(r)\ 2 = N. The solutions of the GP equation 
are obtained numerically by discretizing the wave function ip on a 3D box of size L 
(the number of grid points ranges from 64 3 to 128 3 ). Then, the energy functional 
E[ip] is minimized by using a conjugate gradient algorithm as described in [481 149] . 
for a given realization of the potential and a given particle density n — N/L 3 . The 
numerical solution yields the value of the chemical potential \i as well as the spatial 
particle distribution. Averages over disorder are obtained by repeating the calculation 
for various realizations of the random field. 
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It is useful to rewrite the GP equation (1521) in terms of the dimensionlcss variables 
Udis(r) = VdisW/Vb, f = r/4, V = V /(h 2 /me 2 c ) and Ji = f i/(h 2 /m£ 2 c ). One hnds 
1.=.. 



--V 2 + V o t; dis (r) + ^|0(r)| 2 



(?) = M*) > (35) 



where we have used the wavefunction </>(r) = ^(r) / y^n rescaled in terms of the average 
density n and the healing length £ = 1/ V87ma. From the above equation two regimes 
can be investigated analytically. On the one hand, if the correlation length l c is much 
smaller than the length yj h 2 /mVo fixed by the disorder strength or equivalently for 
weak enough disorder (Vb <C 1), the order parameter is almost uniform |0(f)| 2 ~ 1 
and the effect of the random potential on the chemical potential is just a shift of the 
mean-field energy 

fi = gn + Vo . (36) 

On the other hand, if l c ^> £, one can neglect the kinetic energy term in the GP 
equation (|35[) and use the Thomas-Fermi approximation 

/i = gn(r) + V dis (r) , (37) 

yielding the following result for the local particle density 

n(r) = -[m- Vd i8 (r)]6[/x - 7 dis (r)] . (38) 
5 

Here 0(x) is the Heaviside function: Q(x) = 1 if x > and zero otherwise. The 
normalization condition n = 1/L 3 j drn(r) determines the chemical potential fi in 
terms of the average density n. By using the self-averaging property ([7]) one obtains 
the equation 

£ + = 1 + £ , (39) 

relating the chemical potential to the disorder strength \q. 

i) If Vq <C 5^ then /i = gn + Vo and the disorder acts as a small shift of the pure 
interaction term, similarly to the short-£ c regime of equation (|36[) . 

ii) If Vo ^ gn one finds to the lowest order 

ji = y/2gnV , (40) 
corresponding to an energy per particle E/N = 2[i/3. 

The variation of fi with the disorder correlation length is shown in figure 1121 
where we report the results of the GP equation for the fixed value na 3 = 10~ 6 of the 
gas parameter and disorder strength Vo = 10fcsT°. In this case ^n 1 / 3 ~ 1 and the 
figure clearly shows the two limiting regimes of short and long correlation length. For 
the same disorder strength, in figure 1131 we show instead the density dependence of 
the chemical potential for a fixed value of l c . The equations of state corresponding to 
large and small correlation lengths arc also shown as a reference. 

In the regime t c 3> £ one can make use of the Thomas-Fermi approximation 
for the order parameter which is expected to become more and more accurate as 
£ c increases. Within this approximation one can make contact with the classical 
percolation problem of section [3] by noticing that, if fi is larger than the threshold 
energy e c , the condensate density represented by equation (|38[) is different from zero 
on a percolating path ensuring thus the superfluid behavior of the system. For large 
Vo the chemical potential increases as y/Vo according to equation (|40| . while e c is 
proportional to Vq- The small value of the ratio e c /Vb implies though that the quantum 
phase transition to the insulating state takes place only at extremely large disorder 
intensities. 
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Figure 14. (color online). Equation of state energy vs. temperature of a gas 
with interaction parameter no 3 = 10 -6 and disorder strength Vb = WkgT®. 
Results corresponding to different values of £ c are reported. The limiting cases of 
large correlation length: HF results with LDA (black line) and of small correlation 
length: clean system results shifted by Vo (red symbols and line) arc also shown. 
The open symbols at T = correspond to the results of the GP equation. The 
arrows indicate the superfluid transition temperature. The green line is a T 2 fit 
to the equation of state. 



6.2. Finite temperature 

At T ^ one should combine the GP equation for the condensate with a proper 
description of the thermally excited states in the random potential. The theory 
becomes simple if the disorder is correlated over large distances, as one can apply the 
local density approximation (LDA) within standard mean-field techniques suitable for 
dilute gases. At low temperatures the validity condition requires £ c to be much larger 
than the healing length £, while at higher temperatures the correlation length must 
exceed the thermal wavelength At- 

We use the self-consistent Hartree-Fock (HF) scheme within LDA. This mean- 
field approach is based on the following expression for the elementary excitations of 
the system in terms of their momentum p and position r |50j 



The thermal density of atoms ny(r) is obtained from the momentum integral of the 
distribution of elementary excitations 



The condensate density no(r) is determined by the finite-T extension of the GP 
equation (|38[) in the Thomas-Fermi approximation 



2 



e(P, r) = h Vdis(r) - jti + 2gn(r) . 

Zm 



(41) 




(42) 



™o(r) = - [n - V dis (r) - 2gn T (r)} 
9 



x e[n - V dis {r) - 2gn T (r)} . 



(43) 
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The sum of thermal and condensate density gives the total density n(r), which must 
satisfy the overall normalization 

N = J drn(r) = J dr[n {r) + n T (r)] , (44) 

providing the closure relation of the mean-field equations. The total energy E of the 
system can be calculated from the following integral 

dpdr p 2 /2m f 



{2TThf e c(P,r)/k B T _ 1 

+ | |dr[2n 2 (r)-n 2 (r)]. (45) 

The HF scheme defined by the above equations neglects the quantum depletion 
of the condensate and at T = yields n (r) = n(r), in agreement with 
section 16.11 Furthermore, it neglects the contribution of collective modes (phonons) 
to thermodynamics since all excitations are single particle. This approximation is 
known to be very accurate in dilute non-uniform systems both at high and low 
temperatures !50 . In particular, at low temperatures, one can neglect the thermal 
density contribution to the expression (|41[) of the elementary excitations and one finds 
the simple spectrum e(p,r) = p 2 /2m + |Vdi s (r) — fi(T = 0)|. The term in modulus 
vanishes at the condensate boundaries and thermal excitations accumulate at these 
minima of the effective potential. The single-particle excitations close to these minima 
are the dominant ones at low temperature, being more important than the phonons 
of the bulk condensate. In fact, one can show that at low energy the single-particle 
excitations have the following density of states 

proportional to e 3 / 2 in contrast to g(e) cx e 2 of phononic excitations. A similar 
situation occurs in harmonically trapped condensates [50] . 

The above semiclassical approach provides an estimate of the temperature at 
which Bose-Einstein condensation sets in locally in some deep well. This temperature 
is defined as the point where the local density n(r), corresponding to some deep well 
in the random field, reaches the critical value n(r)A?n = ((3/2) ~ 2.612. By neglecting 
interactions one finds the following implicit equation for the temperature T* 

2vtf \ 3/2 ~ 1 1 



mk B T* J j3/2 1 + jV /k B T* ' (4?) 

The temperature T* is always larger than the temperature T c ° of the occurrence 
of Bose-Einstein in non-interacting clean systems. In particular, for large disorder 
strength one finds T* = T c [C(3/2)V /C(5/2)fc B T c ] 2 / 5 . This effect comes from the 
reduced available volume and the corresponding higher local particle density. In 
the presence of weak interactions, local Bose-Einstein condensation sets in at a 
temperature slightly lower that T*, because density is reduced in the wells of the 
random field due to repulsion and a lower temperature is needed to reach the critical 
value. We would like to stress that the temperature scale T* corresponds to the 
appearance of local condensates at the minima of Vdi s (r) and should not be confused 
with the critical temperature T c at which extended superfluidity sets in. Within 
the above semiclassical approach this latter temperature corresponds to the chemical 
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T=3.2T° T = 2.4T° 




Figure 15. (color online). Particle density profiles in the (x = 0,y,z) plane 
at different temperatures and for a given realization of disorder characterized by 
V = 10fc B T° and t c = O.Gn- 1 / 3 . The T = profile is obtained using the GP 
equation. The average density corresponds to the value na 3 = 10 — 6 of the gas 
parameter. 



potential reaching the percolation threshold of the effective potential Vdi s (r) + 2g r nT(r) 
where, according to equation (|43[) . the condensate density no(r) is different from zero 
on a percolating path. For weakly-interacting systems though, because of the small 
value of the percolating volume fraction $ c , the temperatures T c and T* are very 
close, unless for extremely large disorder intensities. As an example we consider the 
configuration shown in figure HH corresponding to no? = 10~ 6 and Vq = lOfcsTj? in the 
regime of extremely long-range correlation length t c . The value of the temperature 
T* obtained from equation (j47|) is given by T* = 3.63T". The self-consistent solution 
of the HF equations yields a temperature Tbec for the local onset of Bose-Einstein 
condensation in the range 3.5T C ° < Tbec < T*. The transition temperature T c where 
the condensate density percolates is found to be in the range 3.4T C ° < T c < Tbec- 

The HF equation of state and the corresponding transition temperature are shown 
in figure [14] for a fixed density of the gas and for large disorder strength. In the 
figure is also reported the equation of state corresponding to the regime of a very 
short correlation length £ c , where the energy is merely shifted by the average disorder 
intensity Vo from the value of the clean system. This result is consistent with the 
T = prediction (|32|) and the corresponding transition temperature coincides with T c 
in the absence of disorder. 

7. Low temperature thermodynamics 

The energy per particle obtained from PIMC simulations is reported in figure Q3] as a 
function of temperature and for random potentials with different correlation lengths in 
the range of the mean interparticle distance. In particular, for the value l c = 0.67J" 1 / 3 
we also indicate the value of the supcrfluid transition temperature. The equation of 
state and the value of T c corresponding to the limiting regimes of exceedingly long 
and short correlation lengths are also shown as a reference. Three important remarks 
about this figure are worth stressing, a) At fixed particle density and for a given 
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T = 0.8 t£ 




Figure 16. (color online). Particle and condensate density profiles in the 
(x = 0,y,z) plane at two different temperatures, above and below T c . The 
configuration is the same as in figure [T5l 




Figure 17. (color online). Particle and condensate density profiles along the 
(x = 0, y = 0, z) axis corresponding to the configuration of figure [TBI at T = 0.1T°. 



na 3 = 10" 9 na 3 = 10 s na 3 = 10" 3 

Figure 18. (color online). Particle density profiles in the (x = 0, y, z) plane 
for different values of the gas parameter and for a given realization of disorder 
characterized by Vo = 10fc s T° and £ c = 0.6n~ 1/3 . The temperature is T = 0.1T°. 
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Figure 19. (color online). Temperature dependence of the superfluid and 
condensate density without disorder and with strong disorder (Vb = WkgT® and 
ra 1 / 3 ^ = 0.6). The value of the gas parameter is no 3 = 10 — 6 . The grey shaded 
area shows the estimated range of temperatures where the superfluid transition 
takes place. Solid and open symbols refer to two different system sizes to check 
the role of finite-size effects. The solid and dashed lines in the absence of disorder 
are the predictions of a self-consistent mean-field calculation, while the dotted 
lines with disorder are guides to the eye. 




Figure 20. (color online). Dependence on the gas parameter of the superfluid 
and condensate density in the presence of strong disorder (Vq = lOkgT® and 
n, 1 / 3 ^ = 0.6). The temperature is T = 0.1T®. Solid and open symbols refer to 
two different system sizes to check the role of finite-size effects. The error bars 
displayed with the solid symbols corresponding to p s /p come from averages over 
a number of disorder realizations. The dashed lines are guides to the eye. 
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disorder strength the largest suppression of T c is achieved with £ c on the order of 
the intcrparticle distance, b) At low temperature the energetics of the system is well 
described by the GP equation even when the disorder is strong and n 1//3 £ c ~ 1. c) 
For random potentials of large strength and with spatial correlations in the range 
n 1//3 £ c ~ 1, a window of temperatures opens up where the system is normal (i.e. not 
superfiuid), though highly degenerate (nX^ > 1). 

The properties of the "exotic" normal phase displayed at low temperatures for 
strong disorder are worth investigating. We find that the equation of state is well 
fitted by a quadratic temperature dependence as shown in figure Q31 This T 2 -law for 
the energy, and consequently linear dependence of the specific heat, is a remarkable 
feature of the phase. In a clean superfiuid, or in a dirty one with short-range disorder 
correlations (see figure fl"4j). the energy at low temperatures exhibits the T 4 -law typical 
of the thermal excitation of phonons. In the opposite regime of long-range disorder 
correlations one can apply the HF-LDA approach described in section 16.21 which is 
expected to become more and more accurate as £ c increases. Within this approach the 
system consists of large condensate lakes that may or may not be connected through 
a percolating path, the relevant excitations at low temperature arc of single-particle 
nature and are localized at the boundary of the condensate lakes. These excitations 
contribute to the energy with a T 7 / 2 -lavJJ A Bose glass is predicted to exhibit a 
non- vanishing density of states at zero energy [5], which results in a T 2 -law for the low 
temperature equation of state. We interpret the quadratic dependence found in the 
low temperature normal phase as an evidence of the Bose glass phase. 

Another important output of PIMC simulations with the random potential is 
the spatial particle distribution and the distribution of particles contributing to the 
condensate density. In figure [TS] we show the results of the particle density n(r) 
for a given realization of the random potential. Starting from a high temperature 
distribution spread over the entire system, as temperature is decreased, the density 
becomes more and more peaked in correspondence of the minima of the external 
potential. Around T ~ OAT® the system turns superfiuid (see figure ITfjf and the 
density changes only slightly down to the lowest temperatures. Remarkably, the 
comparison with the findings of the GP equation for the same realization of the random 
field is rather good both for the position and the relative intensity of the peaks. 

The comparison between particle density and condensate density profiles is 
reported in figure [16] for two temperatures, above and below T c . The result at the 
lowest temperature shows that the condensate density follows the particle distribution, 
but does not exhibit its pronounced peaks. This behavior is clearly represented in 
figure [171 where we show the profiles along a cut through the plane of figure [16] 
Notice that the particles contributing to the condensate are only a small fraction 
(~ 20%) of the total number of particles (see figure fi"9j). Finally, in figure [TBI we report 
the density profiles as a function of the interaction strength. The figure clearly shows 
the effect of particle localization in the deepest wells of the random field as the value 
of the gas parameter is reduced. 

The behavior of the superfiuid and condensate fraction in the proximity of the 
phase transition with and without disorder are shown in figure 1191 It is interesting to 
notice the large depletion of p s and no even at the lowest temperatures and the fact 
that, in the regime of strong disorder, the superfiuid component becomes significantly 

% In the insulating Bose glass phase the asymptotic low-temperature law is expected to be ~ T 2 
independent of the value of the correlation length. However, for increasing t c , the range of 
temperatures where the asymptotic law applies is suppressed. 
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smaller than the condensate one. Similar results are reported in figure [201 for a fixed 
temperature and for varying values of the interaction strength. The figure clearly 
shows that by making interactions weaker the system will eventually turn normal and 
that the difference between superfluid and condensate fraction disappears with strong 
enough interactions. 

To complete the picture we would like to mention the limit of weak disorder 
(Vom£ 2 /H 2 ) -c 1 and weak interactions na 3 — > considered in Ref. [56]. It explains 
how interactions stabilize the superfluid phase starting from the localized phase of the 
non-interacting gas. At zero temperature the boundary between the insulating (Bose 
glass) and superfluid phases is set (up to logarithmic corrections) by percolation at the 
Larkin energy scale Vq m?£®/h 6 which has to be compared with the chemical potential 
of the weakly interacting gas gn. This leads to the relation for the critical density 
(na) c ~ Vq, i.e. for small Vq only a tiny density of particles is required to off-set 
the localization effects. In practice, one observes extremely robust superfluidity of 
the weakly-interacting Bose gas even in the presence of relatively large disorder, see 
Figs. [JJ HH [20j for na ^> (na) c the critical temperature is essentially the same as in 
the ideal Bose gas. 

8. Conclusions 

We find that in a quantum degenerate bosonic gas a random potential is most efficient 
in suppressing superfluidity if it is correlated over length scales comparable with the 
mean interparticle distance. However, for the typical diluteness conditions of ultracold 
gases, disorder intensities significantly larger than the energy scale fc_eT° set by the 
BEC transition of the ideal gas are required for a significant reduction of the superfluid 
critical temperature and stronger interactions make the superfluid state more robust. 
In the regime of weak interactions and strong disorder, the superfluid transition turns 
out to be well characterized by the existence of a mobility edge, separating localized 
from extended states, that is largely independent of temperature and interaction 
strength. This picture is similar to the percolation threshold of classical particles, 
but we find that quantum localization effects drive the system normal in a large 
region of energies where classical states would be extended. Furthermore, most of the 
particles are localized in the deepest wells of the random potential and only a small 
proportion contributes to the extended condensate state. The effective density of these 
delocalized particles is much smoother due to the screening of the external field from 
the other particles and sets the critical density of superfluidity which however starts 
at significantly lower temperatures compared to clean systems. For larger disorder 
intensities an "exotic" normal phase appears in the degenerate regime, even though 
we can not reach T = 0. This phase is characterized by a peculiar T 2 dependence 
of the equation of state, that is markedly different from the T 4 law of homogeneous 
superfluids and from the T 7 / 2 law found for large condensate lakes within LDA and is 
in agreement with the predictions for the Bose glass phase. Remarkably, some aspects 
of this phase, such as the T = equation of state and the spatial distribution of 
particles can be correctly described using the mean-field GP theory. 
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